The discrete energy method in numerical relativity: Towards long-term stability 
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The energy method can be used to identify well-posed initial boundary value problems for quasi- 
linear, symmetric hyperbolic partial differential equations with maximally dissipative boundary 
conditions. A similar analysis of the discrete system can be used to construct stable finite difference 
equations for these problems at the linear level. In this paper we apply these techniques to some 
test problems commonly used in numerical relativity and observe that while we obtain convergent 
schemes, fast growing modes, or "artificial instabilities," contaminate the solution. We find that 
these growing modes can partially arise from the lack of a Leibnitz rule for discrete derivatives and 
discuss ways to limit this spurious growth. 
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I. INTRODUCTION 



Einstein's theory of general relativity is described by 
a complicated set of coupled, nonlinear partial differ- 
ential equations. Like the Maxwell equations of classi- 
cal electromagnetism, the Einstein equations are over- 
determined, and can be separated into hyperbolic evolu- 
tion and elliptic constraint equations, as in the Arnowitt- 
Deser-Misner (ADM) decomposition [J. The complex- 
ity of these equations is such that analytic solutions have 
been found for only very special configurations, and nu- 
merical studies arguably provide the only means to ex- 
plore a wide varietyof astrophysical and theoretically sig- 
nificant problems [3, y] ■ Unfortunately, these numerical 
solutions are often prone to various instabilities, which 
we divide into three classes: (1) continuum instabilities, 
(2) numerical instabilities and (3) artificial long-term "in- 
stabilities" (ALTI). 

Continuum instabilities exist in the formulation of the 
continuum problem (and thus are naturally reflected in 
any consistent numerical scheme), and are characterized 
by the solution itself or some perturbation of it either 
"blowing up" at a finite time or growing fast in time. 
These instabilities may reflect physical phenomena, such 
as turbulence in fluids or the threshold of black hole for- 
mation [j|, or may be an artifact of the formulation of 
the problem. For instance, ill-posed problems suffer from 
these instabilities, which also manifest themselves as nu- 
merical ones (see 5] , and 6] for an illustration in numer- 
ical relativity). Numerical instabilities are characterized 
by errors that, at a fixed time, increase as the discretiza- 
tion scale is decreased. These instabilities are present in 
many discretization schemes for well-posed problems that 
appear "natural", such as the forward-time, centered- 
space scheme for the advection equation. Finally, ALTI 
are sometimes exhibited even by numerically stable im- 
plementations, whenever errors grow too fast for the de- 
sired simulation time-scale and available computational 
resources. Although these may not be instabilities in the 
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strict sense, given that -as defined here- they are not 
present at the continuum and do go away with resolution, 
they are problematic in practice as they limit the time- 
length of a reliable simulation. In this paper we focus on 
this long-term error growth together with the application 
of the techniques described in '7\ , which throughout this 
work we refer to as Paper I. 

The stability of finite difference approximations 
(FDAs) for simple Initial Value Problems (IVP) with 
constant coefficients and periodic boundaries is often an- 
alyzed in discrete Fourier space, by examining the norm 
of the amplification matrix. This method, however, does 
not generalize in a straightforward way to more compli- 
cated problems, such as those with variable coefficients 
or those with boundaries. A more powerful technique 
is based on the energy method Q, also used to iden- 
tify well-posed continuum problems. This method can 
be used in all problems for which the continuum energy 
method applies, e.g., quasi-linear parabolic and symmet- 
ric hyperbolic partial differential equations with appro- 
priate boundary conditions. The discrete formulation re- 
produces the continuum analysis, including the calcula- 
tion of a discrete energy estimate, which can be used to 
meet sufficient conditions for a stable numerical method. 
Combining these estimates with the consistency of the 
discrete system with the continuum problem, stability of 
the linear problem is equivalent to convergence via the 
Lax equivalence theorem 8]. While our focus in general 
relativity is on the nonlinear Einstein equations, stabil- 
ity of the linearized equations is obviously a necessary 
step for solving non- linear problems, in particular those 
with smooth solutions. In the absence of matter, smooth 
solutions are indeed expected in general relativity when 
appropriate slicing conditions are used PI- 

As described in Paper I, a discrete energy estimate 
can be obtained when using finite difference operators 
that satisfy summation by parts (SEP) [lOj, a discrete 
version of integration by parts, as well as appropriate 
boundary conditions jTj and time integrators [13. Il3l| . 
Paper I presents and extends these techniques in the con- 
text of numerical relativity. In this paper we demonstrate 
these techniques using some common test problems em- 
ployed in the field. We evolve Klein-Gordon scalar fields 



and Maxwell fields on fixed Schwarzschild spacetimes. 
We also consider gauge wave spacetimes, obtained via 
coordinate transformations of flat spacetime. We show 
that we naturally obtain stable and convergent numerical 
schemes by construction, even when an inner boundary 
is used to remove the black hole singularity (black hole 
excision) 1141. 



A second goal of this paper is to understand some 
causes of long-term error growth (ALTI), and ways of 
eliminating or minimizing it by achieving strict stabil- 
ity [ll|. Numerical solutions of the Einstein equations 
are frequently plagued by such long-term errors. Since 
at least some of these appear in numerically stable im- 
plementations, their cause is often attributed to the for- 
mulation of the equations. While this may be the case 
in a number of problems, as discussed below, some nu- 
merical experiments of simple differential equations with 
no continuum instabilities also display this same artificial 
long-term instability — even when discretized in a numer- 
ically stable way — . Hence, it is important to devise nu- 
merical techniques to control this spurious growth in the 
numerical error. We find that this spurious error growth 
is allowed in the absence of sharp energy estimates for the 
discrete problem. This can arise, from among other rea- 
sons, because the discrete derivatives fail to satisfy the 
Leibnitz Rule. With this observation in mind, we dis- 
cuss here a method to remedy some of these problems by 
constructing schemes that suppress excitations of ALTI. 



II. NUMERICAL STABILITY AND THE 
ENERGY METHOD: AN OVERVIEW 

We first review the construction of numerically sta- 
ble finite-difference schemes using the discrete energy 
method for linear, symmetric hyperbolic initial-boundary 
value problems (IBVPs) with maximally dissipative 
boundary conditions. While the discussion here is brief, 
further information can be found in the standard tex t by 
Gustaffsson, Kreiss, and Oliger p\, works by Olsson [ll|, 
and Paper I, for their particular application in numerical 
relativity. SufRcient conditions for a stable discretization 
are fulfilled when: 

1. The continuum problem can be shown to be well 
posed using the energy method. 

2. Spatial difference operators that satisfy SBP on the 
computational domain are used. 

3. Boundary conditions are imposed via orthogonal 
projections pT| . 

4. The semi-discrete equations are integrated with a 
locally stable time i nteg rator, e.g., third or fourth 
order Runge-Kutta [iall3|. 

5. Optionally, explicit dissipation may be added using 
operators that do not spoil the energy estimates, 
and different ways of writing the discrete equations 
may be explored to achieve strict stability. 



This work is organized as follows: We first summa- 
rize, in Section II, some results from Paper I, in order 
to introduce some notation. Section HTll discusses the no- 
tion of strict stability, and we present simple examples 
of "natural" schemes that are numerically stable but not 
strictly stable, for which the errors grow fast in time. 
We then explain the basic strategy to achieve strict sta- 
bility, and show how this eliminates ALTI. In Section Hvl 
we present some test problems from numerical relativ- 
ity: three-dimensional scalar and electromagnetic fields 
propagating on a fixed Schwarzschild black hole geom- 
etry, and one-dimensional perturbations of a dynamical 
slicing of flat spacetime ( "gauge wave" ) . The first cases 
allow us to test black hole excision, while the gauge wave 
has a time dependent background without a sharp en- 
ergy estimate. However, a sharp energy estimate does 
exist for the subsidiary constraint system, and it can 
be used to understand the nature of the ALTI and sup- 
press them. Although more work would need to be done, 
these observations and ideas could be useful when evolv- 
ing Einstein's equations in more complicated situations. 
In Section^we present numerical simulations of the sys- 
tems discussed in Section llVI Finally, in appendix IXI we 
discuss cubical excision boundaries for black hole space- 
times. 



As an example, consider a set of linear, symmet- 
ric hyperbolic equations on a one-dimensional domain, 
X € [o.,b]. This domain is discretized with points Xi = 
a -\- iAx, i — ... TV, where Ax = {b — a)/N, and a 
discrete scalar product. 



N 



{u,v) := Ax y^ cTijUiV 

i,J=0 



3 ' 



(1) 



is introduced for some positive definite matrix with el- 
ements aij. In the limit Ax — > this scalar product 
approaches the continuum one given by 



'dx . 



(2) 



The standard continuum derivative operators and scalar 
products satisfy integration by parts: 



{u',v) + {u,v')^uv\l 



(3) 



Similarly, a discrete difference operator D is said to sat- 
isfy SBP, if there is a scalar product with respect to which 



{Du,v)-\- {u,Dv) =uu|a 



(4) 



holds. The simplest difference operator and scalar prod- 
uct that satisfy SBP on this domain are 



III. ARTIFICIAL LONG-TERM INSTABILITIES 
AND STRICT STABILITY 



D = D+ , (Too = 2 
D = Do, cTu ^1 
D = D^ , ajvN = 



for i = 

for i == 1 ... TV - 1 

for i = A'' 



(5) 



where Dqu = (ui+i — Ui_i)/(2Ax), D-u = {m — 
Ui-i)/ IS.X, D+u = (ui+i — Ui)/Ax, and the scalar prod- 
uct is diagonal: cry = for i ^ j. Higher order opera- 
tors satisfying SBP have been constructed by Strand [13 ■ 
Domains with inner boundaries in several dimensions in- 
troduce additional complexities, and we refer the reader 
to Paper I for more information. 

We sometimes use explicit numerical dissipation. For 
that purpose we modify near boundaries the stan- 
dard Kreiss-Oliger dissipative operator for second order 
derivatives, Qd, so that it does not spoil the semi-discrete 
estimate, i.e., so that it satisfies {u,Qdu) < 0. The op- 
erator Qd, which is fully described in Paper I, is added 
to the right-hand side of the evolution equations. In the 
absence of inner boundaries, it is (on each direction) 



QdUo = -2eAxD\uo, 

QdUi = -eAx{Dl-2D+D^)ui, 

Qdu, = -e{Axf{D+D^fu„ for z = 2, 

QdUN-i = -eAx{Dl-2D+D^)uN-i, 

QdUN = —2eAxD^UN- 



,N-2 



(6) 



We discretize the equations in time using the method of 
lines, and use third-order Runge-Kutta to integrate the 
equations. More precisely, for the ordinary differential 
system 



du 
9t 



L{t, u). 



we integrate the solution m", defined at time i", to the 
advanced time, t""'"-'^ = t" -I- At, obtaining u"+^, through 



7.(1) = 


= «"-K-Atfci, 


«(2) - 


= u" + |Atfc2, 


7."+! = 





3fc2 + 4fc3 



(7) 



where 



fci = i(r, m"), 



L (e + \At, u(iA , 



h = L(r + ^Ai, u(2) 



(8) 



Throughout this paper we call artificial long-term in- 
stabilities (ALTI) those instabilities that are neither 
numerical — the errors at a fixed time do not become 
larger when the resolution is increased, but rather con- 
verge to zero — nor true instabilities of the continuous 
solution. These instabilities, as shown in some examples 
below, can appear quite generically in long-term simula- 
tions. 

There are least two possible ways to suppress ALTFs: 

1. Constructing strictly stable discretizations with re- 
spect to a conveniently defined energy. The equa- 
tions are written such that a sharp continuum es- 
timate also holds for the semi-discrete equations. 
The semi-discrete energy will then closely mirror 
the continuum energy estimate. 

2. Explicit numerical dissipation can be used to elim- 
inate high frequency errors, shifting the spectrum 
of the amplification matrix. 

Although the latter solution is easy to implement, it 
may not be effective if the ALTI excites low frequency 
modes. The first solution, on the other hand, does not in- 
troduce artificial damping, but such a discretization may 
be difficult or even impossible to find for generic systems 
of equations. For example, no general sharp estimate for 
the full Einstein equations is known at this time. Thus, 
these two possibilities can be seen as complementary to 
each other. 

The energy of a solution, w, of a set of linear [i^ partial 
differential equations is defined as 



E{t) = {u,Hu) 



(9) 



where H isa. positive definite symmetrizer for the system. 
This energy |4l| is not unique, reflecting the freedom in 
choosing the symmetrizer, H . The energy method pro- 
vides an estimate for the growth rate of the energy in the 
form of an inequality, say. 



f^-<^) 



(10) 



and thus provides a bound to the energy growth. This 
bound may be much larger than the actual growth rate 
for the continuum solution, in which case we say that the 
energy estimate is coarse. When the estimate is close to 
the actual energy growth rate of a solution, the estimate 
is sharp. While the behavior of the continuum solution 
is independent of a particular bound, as we discuss in 
this paper, numerical solutions sometimes do have a large 
growth if this is allowed. Thus, the sharper the bound is, 
the smaller the possible spurious growth in the numerical 
solution. 

Finding a sharp estimate is problem-dependent, and 
the definition of a convenient energy has to be analyzed 



on each case. One usually attempts to find a sharp en- 
ergy bound by exploiting the freedom in choosing the 
symmetrizer. For example, a well-defined local physi- 
cal energy exists in many physically motivated problems, 
which can be used to select a preferred symmetrizer. The 
evolution equations must be written such that they are 
symmetrizable with respect to this energy, but this is in 
general possible. Although a locally positive definite en- 
ergy for the Einstein equations does not exist, relatively 
sharp estimates exist for some problems |l3, [3 UM ■ 

The energy method can be used to derive discrete es- 
timates. Using difference operators that satisfy SBP and 
representing the boundary conditions through orthogo- 
nal projections, the semi-discrete calculation mirrors the 
continuum analysis. When a discrete energy estimate 
is present, the discrete system is numerically stable by 
construction. A discrete system with the same energy 
estimate as the continuum problem may be defined as 
strictly stable. However, we are particularly interested 
in systems for which the energy estimate is sharp, and 
we therefore reserve the term strictly stable for discrete 
systems that match a sharp continuum energy estimate. 

Assuming a sharp energy estimate exists for the contin- 
uum problem, in this paper we consider how the discrete 
problem might be implemented to preserve this sharp es- 
timate at the numerical level. This can be exploited to 
achieve long-term stability in an otherwise numerically 
stable but long-term unstable scheme. Throughout this 
paper we only use numerically stable schemes, even when 
not explicitly mentioned. Therefore, we concentrate on 
their long-term stability. 

Finally, we note that the sharp energy estimates dis- 
cussed here are for semi-discrete systems, where time is 
continuous but space has been discretized. A sharp esti- 
mate may, in principle, be lost owing to error introduced 
by the temporal integration 53 . The experiments pre- 
sented in this paper and in 20] suggest that (at least in 
the cases considered) if strict stability is lost, the effects 
are not serious. The semi-discrete analysis appears to 
adequately represent the fully discrete simulation. 



A. Example 1: A simple model 

Consider the initial-boundary value problem defined 
by the equation 



u = u' + u/x , 



G[l,2], 



t > 



(11) 



together with some initial data and maximally dissipative 
boundary conditions. Let us consider the numerically 
stable discretization defined by 

ii.^Dui + uJx, , i = O...N, t>0, (12) 

with the discrete derivative, D, given by Eq. Q, third- 
order Runge-Kutta for the time integration [cf. Eq. {T))], 
and maximally dissipative boundary conditions imposed 
through orthogonal projections. Figure ^ shows results 



from simulations using this numerical scheme and u{t = 
0, x) = 1/x as initial data and u(i, 2) = 1/2 as boundary 
conditions, which corresponds to the time- independent 
solution u{t,x) — 1/x . The figure plots the L2 norm of 
the errors with respect to the exact solution as a func- 
tion of time, at different resolutions (obtained on uni- 
form grids, where the number of points ranges from 81 
to 5121). At a fixed time, the errors decrease with reso- 
lution, verifying that, as expected, the scheme is numeri- 
cally stable. At fixed resolution, however, the error grows 
as a function of time, i.e., the solution has long-term error 
growth, or ALTI. 




FIG. 1: This figure shows the L2 norm of the solution error for 
Eq. (1121 as a function of time at many difi'erent discretization 
scales, with initial data it(0, x) = 1/x and boundary condi- 
tions u{t, 2) = 1/2. The solution is calculated with a numeri- 
cally stable discretization that is not strictly stable. Although 
convergent, the error in the solution clearly displays an expo- 
nential growth. The CFL factor used is 0.8. 



To understand this behavior we examine an energy 
norm of the solution. Let the continuum energy | 43| be 
(using H = 1) 



£W^{u,u), 



(13) 



We obtain an energy estimate by taking a time derivative 
in both sides of Eq. (|13|l and using the evolution equation 
(|ll|l . obtaining 

£(1) = (?i,u) + {u,u) 

= {u' + u/x, u) + (u, u' + u/x) 

= {u', u) + {u, u') + {u/x, u) -f {u, u/x) (14) 

= ~u{af + u{hf + {u/x,u) + {u,u/x) (15) 

= —u{af' + - + 2{u/x,u) 

< -- + 2{u/x,u) 

< - + 2a^^{u,u) 



That is, 



~ 4 



(16) 



Only integration by parts was used to derive this esti- 
mate, between Eqs. (|14f) and (|15|l . 

Now consider the discrete problem. We define a semi- 
discrete energy analogous to that one defined in Eq. (|13|l . 
namely 



£;(!) = {u,u) 



(17) 



The derivation of the above continuum energy estimate 
can be repeated at the semidiscrete level, provided the 
difference operator satisfies SBP, since only integration 
by parts was used at the continuum. Therefore, the 
same estimate holds for the semi-discrete energy defined 
in Eq. ^. That is. 



£;(!)< 1 + 2a-ii?(i) . 



(18) 



The continuum estimate, Eq. (|16|) . and the semi- 
discrete counterpart, Eq. H18|l . guarantee that the norm 
of the continuum and numerical solutions will grow with 
a rate bounded by e^*/°. As these are estimates, this 
docs not mean that the exponential growth actually oc- 
curs. Indeed, the continuum solution does not grow, 
since Eq. ()ll|l is the advection equation with the time- 
independent change of variables, v = xu. With this 
transformation, Eq. pi|l becomes v = v' . However, 
the results shown in Fig. ^ indicate that an exponential 
growth does appear in the numerical solution. 

We now define a new energy [ij, that allows for a 
sharper energy estimate and, consequently, a better dis- 
crete system. Let this new energy (defined with H = x^) 
be 



2 



(19) 



As before, an estimate is obtained by differentiating in 
time, substituting in the evolution equations, and using 
integration by parts to get 



^(2) = 



idx = 



x^u 



(u' + -\ dx 



b / ? ? \ ' 
^X U 



dx 



a'^u{a) 



<0 



(20) 



In deriving the estimate H2U|) we have used the Leibnitz 
rule, in the third equality. The analogous semi-discrete 
energy is now 



2 



(21) 



However this semi-discrete energy will not, in general, 
satisfy a discrete version of the continuum estimate given 
in (|20|l . The reason is that although D satisfies SBP, it 



does not generally satisfy the Leibnitz rule, which was 
necessary in deriving the estimate (|2Uf) . It can actually 
be checked that, because of this, discretizing the system 
as in Eq. H12I) will not reproduce the continuum estimate, 
i.e. 



Ei^)^ 



a'^u{ay 



(22) 



An alternative is to rewrite the evolution equation Hll() 
such that the Leibnitz rule is unnecessary to obtain the 
continuum estimate, Eq. (|20(l . Then a semi-discrete ver- 
sion of this estimate does hold if the difference operator 
satisfies SBP. To do so we write Eq. (|llll as 



it = —(ux)' . 

X 

This gives the same estimate as before: 



5(2) 



X uiidx 



xu{xuy dx 



(23) 



\ [a<o)f < 0. 



(24) 
However, in deriving the estimate (|24|l we have only used 
integration by parts, and not the Leibnitz rule. There- 
fore, the semi-discrete energy estimate 



e(^) 



a?'u{aY 



(25) 



follows immediately if the semidiscrete equation is writ- 
ten as 



1 



D{ux). 



(26) 



and D satisfies SBP. 

Figure [21 shows numerical results obtained by evolving 
Eq. (|11|) using two discretizations, given by Eqs. H12I26|I . 
Eqs. Hll|l and H23() are equivalent, they have the same 
set of solutions. However, their discrete counterparts, 
Eqs. (|12|l and (|26|l . do not have the same solutions at 
fixed resolution. In one case the solution error grows as 
a function of time, in the other case it does not. 

It is experimentally found that in this example (and 
others discussed below) the growing numerical error in 
the non-strictly stable discretization is manifest in high 
frequency modes, which can be eliminated by adding a 
small amount of dissipation, as is also shown in Figure El 



B. Example 2: Energy preserving and 
non-preserving system 

Consider now an equation with variable-coefficients in 
the principal part and a lower order term that in principle 
allows for energy growth. 



2cm' 



F{u) 



(27) 



defined on a compact domain x € [a,b]. For simplicity 
we assume that c = c{x) > 0, and define the energy £ = 



10' 



10^ 
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FIG. 2: This figure compares the solution error, at fixed 
resolution, for difi'erent discretizations of Eq. lllll . The dotted 
line shows the one of the lines of Fig.l, i.e., the L2 norm of the 
error for a non-strictly stable discretization, given by Eq. 1121 . 
The solid line shows the error for the same simulation with 
some numerical dissipation (a — 0.01). The dashed line, in 
turn, shows the error for the (non-dissipative) strictly stable 
discretization given by Eq. 1261 . In all these runs Aa; — 1/80 
and CFL=0.8. Even though both the strictly stable and the 
non-strictly but dissipative schemes give long term stability, 
the former discretization is much more accurate. 



(m, u)/2. Homogeneous boundary conditions are given at 
X = b. Rewriting the equation as u = (uc)' + cu' + F(u), 
one can easily show, using only integration by parts, that 
the estimate 

£ = {u\)\l + {F,u) = -uiafcia) + {F,u) < {F,u) 

holds. Therefore, writing the semi-discrete equation as 

ii = D{uc) + cDu + F{u) (28) 

gives the same semi-discrete energy estimate, 

E < (F, u) 

If F = then E,£ < and the continuum and dis- 
crete solutions cannot grow in time. Equations with fea- 
tures similar to those of Ea. H27|) with F — appear in 
the Klein-Gordon and Maxwell fields on fixed black hole 
backgrounds, as discussed later in this paper, where the 
shift plays the role of c. 

If -F 7^ the energy may grow, already at the contin- 
uum; for example, if F = u. In that case one may want to 
discretize in a way such that artificial numerical growth 
is avoided as much as possible. For example, ii F = u, 
discretizing as in Eq. H28|) gives 



E ■ 



-u^{a)c{a)+E 



which is the same equation satisfied by the continuum 
energy. 



Summary 

Fast growing numerical errors can arise even in nu- 
merically stable discretizations. In some cases they can 
be traced to the failure of the Leibnitz rule for discrete 
derivatives, and failure of sharp continuum energy es- 
timates to hold in the discrete problem. The use of 
difference operators satisfying SBP, appropriate numer- 
ical boundary conditions, and rearrangement of terms, 
give raise to strictly stable schemes which satisfy sharp 
semidiscrete energy estimates, if such estimates are avail- 
able at the continuum. 

The goal of these rearrangement of terms is actually 
not only to avoid the Leibnitz rule. Clearly, since the 
growth rate for the energy usually involves derivatives, 
in general there would still be a discrepancy between the 
analytical and numerical estimates even if the Leibnitz 
rule was satisfied, arising from the difference between 
numerical and analytical derivatives. These differences, 
however, can also be controlled by appropriately rear- 
ranging terms (see, for example. Sections II V A III V A 21 
andRefs. SEllll). 

In the following we discuss more complex examples: 
the propagation of scalar and electromagnetic fields on 
the Schwarzschild spacetime in three dimensions, and 
one-dimensional linear perturbations of gauge waves. In 
the former cases, we do have readily available sharp en- 
ergy estimates derived from physical principles. In the 
latter a sharp estimate is available for the subsidiary sys- 
tem describing the constraint evolution. 



IV. FORMULATION OF SAMPLE PROBLEMS 

In this section we describe some model problems com- 
monly used in relativity, while corresponding numerical 
results are presented in the following section. The first 
two problems are Klein-Gordon and Maxwell fields on a 
Schwarzschild black hole background, and linearization 
perturbations of the gauge wave is the third one. 



A. Scalar and Maxwell fields on a black hole 
background. 

Black holes contain physical curvature singularities 
that can not be included on the computational domain. 
Black hole excision eliminates this problem by remov- 
ing a region around the singularity from the computa- 
tional domain [lj| . This region must lie within the event 
horizon, and be causally disconnected from the spacetime 
outside the black hole so that the excised region does not 
affect the physics outside the hole. As discussed in the 
appendix, this region must be carefully chosen for these 
conditions to hold. 

Numerical implementations of excision have tradition- 
ally presented a challenge in deciding how points near the 



inner boundary are updated. Especially in higher dimen- 
sions, many combinations of one-sided derivative stencils, 
interpolation, and extrapolation have been investigated, 
mostly through experimentation [2II2I I2I I2ll27ll23| 
(for an alternative excision strategy using pseudospec- 
tral methods see |23)- However, as explained in paper 
I, derivative operators that satisfy SEP together with a 
diagonal inner product uniquely specify the excision al- 
gorithm. The procedure is unambiguous, and for linear 
problems guarantees a stable implementation. We follow 
this procedure, detailed in Paper I, to excise the singu- 
larity. 

We consider three-dimensional evolutions on a 
Schwarzschild background, in either Painleve-Gullstrand 
(PG) and Kerr-Schild (KS) coordinates. The four- 
dimensional metric, written in terms of ADM quantities. 



ds^ = -a^dt^ + /ly (dec* + P'dt){dx^ + P^dt) . (29) 
In PG coordinates 



a = 1 , (3' ^ 



2m\ 



1/2 



This formulation has a number of desirable features 
when the background metric is stationary J2j|. For ex- 
ample, Eqs. l|5^ are symmetric with respect to the phys- 
ical energy, and a straightforward estimate shows that 
this energy does not grow when maximally dissipative 
boundary conditions, which are automatically constraint- 
preserving for this formulation, are given. Moreover, 
Eqs. H32|) have been written such that the Leibnitz rule 
is not needed to guarantee that the energy will not grow. 
For this reason, we refer to the discretization obtained 
by replacing di by Di in Eqs. (|32|) as energy preserving 
or strictly stable. The physical energy, however, is not 
positive definite inside the black hole, since the Killing 
vector becomes space-like there. We nevertheless choose 
this energy-preserving discretization, since it guarantees 
long-term stability if the computational domain does not 
include the black hole. More sophisticated formulations 
and discretizations, which are globally symmetric hyper- 
bolic and conserve the physical energy in the exterior 
region, up to a distance arbitrarily close to the event 
horizon, can also be constructed 



,riiy 



f^ij — Oij , 



(30) 



2. Maxwell fields on a Schwarzschild background 



while in KS coordinates 



a 



r 



f3' = 



r + 2M 
2M , 



1/2 



::X 



h-ij — dij 



2M 
2M x'x' 



(31) 



where r = {x"^ + y"^ + z^)^!"^ and M is the black hole's 
mass. We excise a cubical excision region centered on 
the black hole, and require the length of the cube, L, to 
satisfy L < 4M/(3v/3), or L < 0.78M. Larger cubes, as 
discussed in the appendix, result in excision boundaries 
that are not purely outflow. 



The Klein-Gordon field on a Schwarzschild background 



The massless Klein-Gordon scalar wave equation, 
3°''VaVb<i> — 0, where g'^^ is the inverse four-dimensional 
metric, can be be written as 



5t$ 
dtU 

dtdi 



n, 

+h-^/^ddah^'^H'^dj) 

9,n, (32) 



where we have introduced new variables 11 := 9t$ and 
di := di^ to write the system in first order form. Addi- 
tionally, we define H^^ := /i*-' - 
inverse of the three-metric, hi 



a P^(3^ , where /i*-' is the 
, and h = det{hij). 



As a second example of classical fields on a black hole 
background we choose the vacuum Maxwell equations. 
We write the Maxwell equations in terms of the Faraday 
tensor. Fab, and its dual 



Vr„,K 



be] 



Vr/F, 



he] 



(33) 



The dual of the Faraday tensor is *Fab '■— hsab'^'^F^ 



cd. 



where Sabcd is the completely anti-symmetric symbol. 
Contracting the above equations with a time-like vector 
field u° we get the evolution equations: 



f^u^ ab 
*^u ^ ab 



-2V[,(Ffc],u^), 

-2\/[aCFb]cU^). 



(34) 



Defining a foliation by dragging an initial space-like 
hypersurface along the integral curves of the vector u°- 
(each labeled with the proper time of the integral lines 
of u", t), and pulling back the above equations onto each 
of these hypersurfaces we get 



Cu(i>i,{Fab) = -2d[a(t>*{Fb\cU''), 
CuM*Fab) = -2dlaM*Fb]cU'), 



(35) 



where (f>* is the pull-back map into the hypersurface. No- 
tice that only exterior derivatives appear in the equa- 
tions, so the geometry only enters through the star op- 
eration. The electric and magnetic fields are typically 



defined as E„ 



Fabn , and Ba 



bed 



FbeUd, allow- 



ing one to decompose Fab and *Fab into {n'^ defined by 

Eab'' BcUd, 



Fab 



-2E\ 



anb] 



*Fab = 2B[anb] + Eab EcUd 



(36) 



Contracting Eqs. 135() with ^e^'^'^'^Ud, we get 

C^E, = s,^''djiaBk-Ski„,P''E'^) (37) 

CuB, = ~e,'''d,{aEk+ekirnl3'B"'). (38) 

As with the Klein-Gordon equation, this system is sym- 
metric hyperbohc outside the event horizon, where u" is 
time-hke, and only strongly hyperbolic inside the event 
horizon, where m° is space-like, since the symmetrizer is 
not positive-definite there. 

We supplement the above equations with a no- 
incoming radiation boundary condition. No condition 
is needed to preserve the constraints, for they propagate 
along M°, as are the outer boundaries. Thus, we set the 
incoming modes to zero at a boundary with normal rua, 
via the projector defined as 



PiE,B) 



= I-Y,e,e\E,B) 

A+ 



— {Ei , B, 
1 



{E, - e^'^B^mk, B, + e,^'=^,mfc)(39) 



where e^ is the eigenvector basis, 
define 



is the co-basis, we 



E, := E,~{m,-l3^)m^Ej/{p^mi + l), 
B, := B^~{m^-f3^)m^BJ/{p^ml + l), 

and where the sum extends only over eigenvectors with 
positive eigenvalues. 

As with the scalar field, we ask whether the Maxwell 
equations can be written such that the discrete problem 
has a sharp energy estimate, giving numerical solutions 
without ALTI. The background spacetime is invariant 
under time translations, implying a conserved quantity 
or energy. This energy is the integral of Tabn°'u\ where 
w° is the killing vector field (timelike outside the black 
hole) over the hypersurface. For the PG foliation this 
energy is 



^K = \ f[E'E, + B'B, + 2E'B^p''e^,k] dV. 



(40) 



Just as with the Klein-Gordon equation, this expression 
is not positive-definite inside the horizon, refiecting a 
property of the geometry, not of the equation or coor- 
dinates used. 

To discretize the above system we use the method of 
lines and substitute in Eqs. (|37I38|I the partial deriva- 
tives with difference operators satisfying SBP. In this way 
strict stability with respect to the discretized version of 
the above energy is attained. 



B. Linear gauge wave propagation 

One simple and common test for numerical implemen- 
tations of the Einstein equations is a gauge wave defined 



by 

which describes flat spacetime with a sinusoidal coordi- 
nate dependence, of amplitude A, along the x direction. 
The non-trivial ADM variables for this metric are 



!Jxx 



^A sin(7r(a: — i)) 

A 



(42) 



Kxx = -^ cos (^(x-i))e-4/2 -"«--*», (43) 



together with the gauge condition 

^ ^ ^A/2sin{ii{x-t)) 

P' = 0. 



(44) 
(45) 



We evolve the Einstein equations u sing the symmetric 
hyperbolic formulation presented in |3l| with dynamical 
lapse and the time-harmonic gauge. The formulation is 
cast in first order form by introducing the variables Ax '■= 
dxCt/a and dxxx '■= dxgxx- Rather than solving the full 
Einstein equations here, we study linear perturbations of 
a background metric given by Eq. H41() . Furthermore, we 
simplify the system by allowing perturbations depending 
only on (i,x), in the form 

9xx — 9xx ' OQxxi 
^xx — ^xx ~r OPi^xxi 

^XXX ^XXX I ^^XXX7 

a = a + 5a, 

•^X -^X I ^-^x ■ 

This restriction can be justified by noting that full non- 
linear evolutions of data defined by Eq. H41|l show that 
only these variables vary dynamically [2C| , and variations 
in all other variables are consistent with round-off error. 
The resulting linearized equations are (henceforth drop- 
ping the S notation) 

Att 
dta = -ATTC0s{4>)a- Kxx + -;r:rC0s{4>)gxx,{^^) 

la 

1 All 

dtAx = ~—dxKxx + TTT cos{4>)Kxx , 





a la 






A-k"^ 
~'2^ (Acos(0)^ + sin(0)) gxx 






+ -^coa((j))dxxx 






— 7rcos{(f>)Ax + -rp- sm 0)a , 
I la 


(47) 


dtQxx = 


= —AaiT cos{(f>)a — 2aKxx , 


(48) 


dtKxx -- 


= -adxAx — cos(0) A 

— (-2sin((/)) + Acos(0)^) a 

-~ATrcos{(j))K^^ 






+ -rrCos((j))dxxx , 
4a 


(49) 



Ot^xxx 



AaiT^ (sin((/)) — ^cos(0)^) a 

—Aa-KCOs{(l))Kxx — A6?T:cos{(f))Ax 
-2adxKxx , (50) 



where we have defined (j) '■— '^{^ ~ t)- This system is 
symmetric hyperbolic; three characteristic speeds are 
and the other two ±1. 

We consider here an initial boundary value problem 
on the domain x E [—1/2,3/2]. We implement bound- 
ary conditions via the orthogonal projector method de- 
scribed in appendix D of Paper I. The positive/negative 
eigenmode decomposition of the matrix HA is W± = 
— {Ax^Kxx) (see Paper I). Here, the matrix A is the prin- 
cipal part, and H is its symmetrizer. For one-dimensional 
problems, H = {P^^)'^P^^, and P is the matrix that di- 
agonalizes A, i.e., P^^AP = A, with A a diagonal matrix. 

The system must satisfy two constraint equations, cor- 
responding to the definition of the variables dxxx and Ax ■ 
When linearized, these constraints are 



= Cx 



-dxdx 



(51) 



1 / Att \ 

= Ca=Ax~^ idxa - — cos(0)a j . (52) 

The other constraints are automatically satisfied owing to 
the restricted form of the allowed perturbations. Despite 
the simplicity of this system, the numerical integration of 
these equations reveals some interesting features. Indeed, 
straightforward discretizations of this system have lead 
to exponentially growing solutions with quickly growing 
constraint violations. That this instability has been ob- 
served by different groups in a variety of situations, led 
to the speculation that constraint violations drive the in- 
stability. We offer here a different interpretation of the 
results by examining the evolution of the constraints. 

Evolution equations for the constraints can be obtained 
by taking the time derivative of Cx and C^, and substi- 
tuting in the Einstein equations, giving 



dtCx — —A-Kc? cos(0)Ca , 



(53) 



JK'K A.'K 

dtCA = ^ cos(0)C, - — cob{4>)Ca ■ (54) 

These equations can be integrated in closed form to yield 

Cx = G{x)a^ + H{x)a , (55) 



Ca = G{x) + 



H{x) 
2d ' 



(56) 



where G{x) and H{x) are free functions determined by 
the initial data: 

G{x)a' = 2a^CA{<^.x)-Gx{Q,x), 
H{x)a = 2{Cx{0,x)-a^CA{0,x)) . 

Clearly the constraint behavior is purely oscillatory. 
Hence, any growth in the constraints observed during 
evolution must be purely numerical. We will show in the 
next section that this is indeed the case, and that the 
error growth can be eliminated. 



V. NUMERICAL SIMULATIONS 

We now present some numerical simulations of the 
model problems discussed in the previous section. For 
the black hole tests, we examine configurations where 
the computational domain is outside the black hole and 
others that do contain an excised black hole. The former 
configuration allows us to investigate the evolution of a 
symmetric hyperbolic system with boundary interactions 
without the difficulties introduced by excision. We then 
examine the more demanding case in which excision is 
required, and the evolution systems here considered are 
only strongly hyperbolic inside the black hole. 



A. Wave propagation on a Schw^arzschild 
background in KS coordinates 

This section presents numerical results for the massless 
Klein-Gordon scalar field on a Schwarzschild black hole 
background. We present some tests of our code, and 
then discuss differences between the strictly stable and 
non-strictly stable discretizations. 



1. Case 1: Computational domain outside the black hole 

Consider an uniform computational domain that does 
noi contain the black hole, {x,y,z) e [1.5M, S.SAf]'^. Ini- 
tial data with compact support are given to H, the time 
derivative of the scalar field, as 

n = 10^[(a;-3Af)(a;-4M)(y-3M)(y-4Af) 

X {z - 3M){z - 4:M)f (57) 

when {x,y,z) e [3M, 4Af]'^, and zero otherwise. All 
other fields are initially set to zero [i^. This data is 
then evolved using two possible discretizations. The first 
one is that one obtained by direct substitution of partial 
derivatives by discrete ones in Eq. 1)32(1 , this is the strictly 
stable discretization. The second one is obtained by ex- 
panding all derivatives of products with the product rule 
and then replacing partial derivatives of field variables 
by their discrete counterparts. We refer to this as the 
"naive" discretization. 

Both discretizations were used to evolve for around 250 
crossing times, under different choices of Courant factor 
(A = 1/4, 1/2, 1). In all cases we found no sign of long 
term growth when monitoring the variables, though a 
closer inspection reveals several salient features of the 
obtained solutions. Namely, a self-convergence analysis 
show results consistent with second order accuracy for a 
few crossing times (with better values for the smallest A 
). However, at later times, the observed convergence rate 
deteriorates considerably and achieve artificially large 
values. The time at which this behavior starts coincides 
with the fields mostly leaving the computational domain. 



10 



For instance, Figure |31 displays the L2 norm of 11 ob- 
tained with the strictly stable and "naive" discretiza- 
tions. The figure shows the L2 norms of 11 calculated 
at two resolutions, for 250 crossing times. The large 
decrease in the norm is a sign that basically most of 
the fields leave the domain after a few crossing times. 
Both discretizations give a long-term stable algorithm 
for this case. But, as we will discuss later, the naive dis- 
cretization without dissipation becomes long-term unsta- 
ble when the black hole is included on the domain. 
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FIG. 3: This figure shows the L2 norm of 11 in a scalar 
field evolution outside the black hole for about 250 crossing 
times. The figure compares results for a naive discretiza- 
tion (see text), and a strictly stable discretization at two 
resolutions, A = A//10 and A = M/20. The domain is 
{x,y,z) £ [1.5,5.5], A — 0.5, and no dissipation is added to 
the solution. 



and both discretizations (strictly and non-strictly stable) 
are evolved. Figure ^ shows the self-convergence factor 
(as before, in the I/2 norm) of $ and 11 versus time for the 
strictly stable discretization (similar factors are obtained 
for the non-strictly stable one) at the resolutions A — 
M/10, Af/20 and M/40. Convergence factors close to the 
expected value of two are observed, with no qualitative 
difference in the convergence factor as time progresses. 
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FIG. 4: This figure shows the self-convergence factor as a 
function of time for the fields 11 (solid line) and $ (dashed 
line) on a domain outside the black hole, using the strictly 
stable scheme (similar factors are obtained for the naive one), 
for about 25 crossing times. These factors are calculated using 
resolutions of M/10, M/20, and M/40. The obtained factors 
are close to the expected value of two. The Courant factor is 
A = 1, and no dissipation is used. 



To verify that we get the expected self-convergence fac- 
tor in long term simulations when the fields do not decay 
to very small values after having left the domain, we now 
show results of a complementary set of tests on the do- 
main X e [2.1M,6.1M], (j/,z) e [-2M,2M]. Here we 
concentrate in what corresponded to the worst case sce- 
nario considered, given by A = 1. We set the incoming 
modes on all outer faces to zero, except on the x = 6.1M 
face, where the time derivative of the incoming fields is 
set to 



VF+ = < 



(y2 + z2 - Mf{M/t - If sin(i) 

ifi> l,(y,z)e [~M,M] 
otherwise 



This boundary condition is imposed through an orthog- 
onal projection, as explained in Paper I. Initial data is 
given by 

n = IQ' [{x - iM){x - AM){y - M){y + M) 

X {z - M){z + M)]'^ , (58) 



2. Case 2: Domain with excision 

In this section, the computational domain contains the 
black hole, which is excised with an inner cubic boundary 
centered on the black hole. The domain is defined on 
[— 4.5M, 4.5M]'^, and the excision cube has a total length 
of 0.6M. Thus, the faces of the inner boundary are at 
±0.3M. 

We first verify that the excision algorithm does not dis- 
play ALTI when employing the strictly stable scheme. As 
mentioned earlier, very little freedom is available in con- 
structing difference operators that satisfy SBP. Indeed, 
for our second-order code with a diagonal scalar prod- 
uct, the stencils arc uniquely specified. Figure [3 shows 
the L2 norms of 11 and $ for a long -up to 10'*M- evo- 
lution (with no dissipation), designed to detect slowly 
growing modes. As shown in this figure, the algorithm 
specified by requiring SBP and the equations arranged so 
as to preserve the energy is long-term stable even with 
excision boundaries. For this run non-trivial initial data 



11 



are given, as before, only to H. 



j-8^-1 W^ _ ^ ^2 _ s^y (rO - r) if \r - ml < \S\ 
11= ( (59) 

otherwise . 



with rg = 2.5M and 6 = 1.5M. When studying the con- 
vergence of the code one again faces the fact that after a 
few crossing times the convergence value become mean- 
ingless. The behavior observed is analogous to the case 
where the black hole was outside the computational do- 
main. Namely, close to second order convergence is seen 
during the first few crossing times, and unrealistically 
high values are obtained afterwards. In Fig. we show 
a short-term convergence test (about two crossing times) 
with excision. Fig. [7| shows the corresponding L2 norms 
for $ and H. 



Convergence factors, wave equation with excision 



Long term excision run 
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FIG. 5: This figure demonstrates the long-term stability of 
the strictly stable discretization of the scalar wave equation 
when an excised black hole is on the grid. The I/2 norms 
of n (solid line) and $ (dashed line) are shown for about 
2500 crossing times. No slowly growing unstable modes are 
detected in these runs. This run was performed with 81^ grid 
points, A = 0.25, and no dissipation is used. 



As in the previous case, we also compare results ob- 
tained with the strictly stable and the naive discretiza- 
tions. Long term stability was found for both cases when 
the computational domain was outside the black hole. 
With excision, however, the solution obtained with the 
naive discretization quickly becomes corrupted, as shown 
in Fig.|Sl As shown in this figure, the amplitude decreases 
with resolution, indicating that the growth is spurious 
and that the naive discretization suffers from ALTI. The 
strictly stable discretization, however, with a sharp semi- 
discrete energy estimate, remains well-behaved. The 
ALTI that contaminate the naively discretized solution 
are found to be of high-frequency type, which can be 
managed by adding explicit dissipation. Figure El shows 
results using the same initial configurations, except that 
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FIG. 6: This figure shows a self-convergence test for 11 (solid 
line) and <E> (dashed line) using the same configuration as in 
Fig. |5] The self-convergence test is performed on uniform 
grids with 81^,161^, and 321'' grid-points. The convergence 
factor is close to the expected value of two. The measured con- 
vergence factor drops slightly from near two after one crossing 
time, when most of the field has left the domain and after- 
wards increases to unrealistic high values. 



now dissipation is added to the naive discretization. Now 
all of the runs produce long term stable simulations. 



B. Maxwell fields in a Sch'warzschild background in 
PG coordinates 



As with the scalar field tests above, we perform tests 
with the Maxwell equations on two kind of domains. 
Specifically, we choose 

1. Computational domain outside the black hole: 

y,ze [-4M,4M]2 and x G [3M, IIM]; 

2. Computational domain including part of the black 
hole: {x,y,z) € [-8.5M,8.5Mf. The singularity 
is excised by means of a cubical region defined by 
ix,y,z) e [-0.35Af,0.35M]3. 

The initial data describe a "pulse" of electromagnetic 
field given by 



F,.= 



yS-^ [{r - ro)2 - S^]\,,s{x^ - xi) if |r - ro| < <5 
) otherwise , 



where Fi stands for Ei and Bi, a is a strength parameter, 
Xq is the location of the pulse's center, and 6 its width. 
We discretize in a strictly stable way, as explained in 
Section irvT2l 
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FIG. 7: This figure shows L2 norms for the runs shown in 
Figini The relative difference in the norms for <1> at any two 
different resolutions, is at most about one percent. Similarly 
for n, except at the last point shown, where this field has 
mostly left the grid. 



1. Case 1 (domain outside the black hole) 

Let a = 1, Xq = (7M, 0, 0) and 5 = M, and we set the 
dissipation to zero. Figure ITUl shows ||i?2:||2 vs. time for 
four different resolutions (81^,161^,241^ and 321^ grid- 
points, corresponding to A = M/10, M/20,M/30 and 
M/40, respectively), while the inset shows the decay in 
time for a long-term run with the coarsest resolution. 
Figure ^J shows the self-convergence factor of E^ for 
a short time, computed with 81'^, 161^ and 321'^ grid- 
points, and 161^, 241^, and 321^^ grid-points. As ex- 
pected, the factor gets closer to two when resolution is 
increased. 
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FIG. 8: This figure shows the L2 norm of 11 in a scalar field 
evolution on a domain containing an excised Schwarzschild 
black hole. The solid and dot-dashed lines show the L2 norm 
of n obtained using the naive discretization of the wave equa- 
tion. Although the amplitude decrease with resolution, this 
solution quickly becomes long-term unstable, indicating an 
artificial instability. The strictly stable discretization (shown 
with dashed and dotted lines), however, remains long term 
stable and well-behaved. The results are calculated with the 
resolutions A = M/9 and A = A//18, A = 0.5, and no dissi- 
pation is added. 



times, even with a strictly stable discretization. The rea- 
son for this appears to be that the energy inside the black 
hole is not positive definite and boundedness (either at 
the continuum or numerical level) is not guaranteed (it 
is not clear why this divergence does not appear in the 
scalar field case, when the singularity is excised, since the 
energy is also non-positive definite in that case.). High 
frequencies are seen in the solution, and we therefore add 
a small amount of dissipation, a = 0.05. Figures [T^ and 
1131 show the norms and convergence data similar to the 
figures for the domain outside the black hole. As be- 
fore, the simulations are long-term stable for long times. 
The self-convergence factors for the coarsest resolutions 
drop quickly from their initial value, indicating that the 
coarsest solutions are not in the convergent regime. In 
addition to effects caused by the fields leaving the compu- 
tational domain, using a cubical excision boundary forces 
us to excise very close to the black hole, where the space- 
time and field gradients are very large, effectively requir- 
ing higher resolution. 



C. Linearization around the gauge wave 



2. Case 2: Domain with excised black hole 

Now let a = I, xl ^ (6.5Af, 0, 0) and 5 = M. In con- 
trast to the scalar field runs, solutions of the Maxwell 
equations without explicit dissipation diverge at long 



The gauge wave is a useful test problem for examin- 
ing some of the difficulties found in the simulation of the 
Einstein equations (see [a, IH, 34, 36] ), and to illus- 
trate the usefulness of the techniques presented in Pa- 
per I. To this end, we present results obtained using a 
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FIG. 9: This figure shows the L2 norm of 11 in a scalar field 
evolution on a domain containing an excised Schwarzschild 
black hole. The solid and dot-dashed lines show the solu- 
tions obtained using the naive discretization with dissipa- 
tion, the latter eliminating the artificial instability noted in 
Fig. I2I The strictly stable discretization (shown with dashed 
and dotted lines), without dissipation, is shown for compari- 
son. The results are calculated with the resolutions A — M/9 
and A = M/18, A = 0.5, and the dissipation parameter is 
a = 0.02. 
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FIG. 10; This figure shows the L2 norm of Ex vs. time in an 
evolution of the Maxwell fields at four different resolutions for 
a strictly stable discretization. The computational domain is 
outside the black hole. The inset shows a long-term evolu- 
tion, 3000M, using the coarsest resolution. The runs were 
performed on uniform grids with resolutions Af/10, Af/20, 
M/30 and A//40, and the Courant factor is A = 0.25. 



FIG. 11: This figure shows self-convergence factors for Ex for 
the runs of Fig. 1101 The factors improve with resolution, in 
the sense that they closer to the expected value of two. The 
labels 81-161-321 and 161-241-321 indicate the number of grid 
points in each direction used to perform the self-convergence 
test. 
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FIG. 12: This figure shows the L2 norm of Ex vs. time in 
an evolution of the Maxwell fields on a computational do- 
main containing an excised black hole for a strictly stable 
discretization. The norms are shown for four different resolu- 
tions, and the inset shows a long-term evolution, 3QQQA'f , at 
the coarsest resolution. The runs were performed on uniform 
grids with resolutions M/10, M/20, M/30 and A//40. The 
Courant factor is A = 0.25, and dissipation is added to the 
evolution equations, with a = 0.05. 



straight-forward discretization of the right hand side of 
Eqs. (|lS|) - (|5n|l . and show results obtained for long term 
evolutions with and without dissipation. Moreover, we 
consider cases both where the constraints are initially 
satisfied, and where they are not. The latter illustrates 
that even when large constraint violations are introduced 



in the initial data, the solutions to the linearized equa- 
tions remain well behaved. This, of course, may not nec- 
essarily be the case for the non-linear Einstein equations, 
as for these relatively large amplitudes the linearized sys- 
tem need not be a close approximation to the non-linear 
ones. Or even if the amplitudes were small, there could 
be purely non-linear effects in the solution which would 
naturally be missed in a linearized analysis |35| . 
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FIG. 13: This figure shows convergence data for the runs of 
Fig. 1121 The self-convergence factors are calculated from Ex- 
The excision boundary is very close to the black hole, and the 
large gradients encountered in the fields and spacetime near 
the excision boundary require very high resolution to resolve 
them adequately. The coarsest solutions are not resolved well 
enough to be in the convergence regime. However, it can 
be seen from the Figure that the convergence improves with 
resolution. The labels 81-161-321 and 161-241-321 indicate 
the number of grid points in each direction used to perform 
the self-convergence test. 



for longer times with increasing resolution, this approach 
is impractical for long runs. Figure [T^ shows the conver- 
gence (to zero) factors in the L2 norm obtained for Cj,- 
For the highest resolutions, the factor is very close to its 
expected value. 

An exponential behavior for gxx is also seen in Figure 
1161 Even though in principle this quantity could already 
grow at the continuum, the fact that at fixed time the 
norm of g^x decreases with increasing resolution suggests 
that, as discussed below, this is an artifact of an ALTI. 
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For these tests we adopt trivial initial data for all fields, 
except for gxx, dxxx and Kxx, which have compact sup- 
port in a; € [0, 1] and are given by 



FIG. 14: Logarithm of the L2 norm of Ca vs. number of 
crossing times for four different resolutions. After an initial 
stage where the constraints are well behaved a clear exponen- 
tial mode dominates the solution. 



9xx 


— Ksin(x) (a; - 


-0.5)2- 


-0.25]\ 


(60) 


^xxx 


= fJ-dxgxx, 






(61) 


Kxx 


= --i^Oxgxx, 






(62) 



where k, /i, and ^ are constants. We set A = 0.1 in 
the background solution and the computational domain 
is taken to be a; G [—0.5,1.5]. Boundary conditions are 
applied via Olsson's projection operators. Here, as we 
are interested in the long term behavior of the solution, 
we couple the incoming modes the outgoing ones such 
that 

W+ = 0.5W. . 

For cases where the constraints are initially satisfied, we 
choose /i = 1 = ^, K = 0.01; otherwise we set /i = 0, k = 
0.01 and ^ = 1. In all these simulations a Courant factor 
of A = 0.25 is used, to compare with similar simulations 

UMM M 

Figures ri4II6l show results for constraint-satisfying ini- 
tial data, without dissipation. Figure ^^ shows the L2 
norm of the constraint, Cj[ and its convergence to zero 
with increasing resolution. Initially the norm Ca is al- 
most constant. However, at later times a clear exponen- 
tial mode is observed. While the solution is well-behaved 
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FIG. 15: Convergence factors (to zero) for Ca as a function 
of number of crossing times, comparing runs with 81 and 161 
points (dotted line), 161 and 321 points (dashed line) and 
321 and 641 points (solid line). Notice that even when the 
solutions are growing exponentially (see Figs ll4ll^ . the con- 
straints converge to zero. 
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FIG. 16: Logarithm of the L2 norm of g^x vs. number of 
crossing times for four different resolutions with no dissipa- 
tion. At early times, the norm remains fairly constant but at 
later ones an exponential growth dominates the solution. 
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FIG. 18; 1/2 norm of Ca vs. number of crossing times for 
four different resolutions. Convergence to zero is observed 
and only slow (linear) growth is observed. 



The ALTI for C_4 are, as in the previous examples, 
found to be of high frequency type and, therefore, the 
addition of dissipation eliminates these growing modes. 
Figures [T7I 1201 iUustrate the behavior of simulations with 
the same initial data and boundary conditions as those 
of Figures I14I16I , except that now some dissipation is 
added, a — 0.01. The norm of Ca is illustrated in Fig. 1181 
where now only a slow growth is observed at late times. 
Fig. 1201 shows, as before, the convergence factors (to 
zero) for C^, in the Li norm. Finally, figure [T7I shows 
the hi norm of g^x for different resolutions. Only a linear 
growth is observed in the solution is found; furthermore, 
this growth decreases with resolution. 
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FIG. 19: Convergence factors calculated with Ca vs. num- 
ber of crossing times comparing runs with 81 and 161 points 
(dotted line), 161 and 321 points (solid line). 



As a last example we consider initial data violating 
the constraints. Here, similarly good behavior is ob- 
served when dissipation is employed, in the sense that, 
as predicted by the continuum analysis of Section IIVBI 
the constraints oscillate in time. For instance, figure 1201 
shows the Li norm of the constraint Cx vs. time. As 
expected, these converge to a roughly constant behavior 
as resolution is increased, and no growth is observed. 



FIG. 17: 1/2 norm of Qxx vs. number of crossing times for four 
different resolutions. Only a small growth is observed in the 
solution. 



VI. CONCLUSIONS 

The construction of stable numerical methods for IB- 
VPs often involves a mixture of intuition and experimen- 
tation. The choice of numerical scheme and boundary 
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FIG. 20: 1/2 norm of Cx vs. number of crossing times for 
three different resolutions, corresponding to the case where 
the initial data does not satisfy Cx ~ 0. As expected from 
the continuum analysis, even in this case the norm does not 
grow in time. 



in handling otherwise artificial long term instabilities. 

It is expected that these techniques will be important 
in a number of fronts, in particular in the simulation of 
Einstein's equations. 
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conditions are often considered as necessarily separate 
ingredients. For certain problems, however, the energy 
method for discrete systems is an integrated, systematic 
method for creating numerically stable schemes and an- 
alyzing their behavior. While this method has not yet 
been widely used in general relativity, many current out- 
standing problems in the field vitally depend on under- 
standing the subtle interplay of boundaries with the inte- 
rior computational domain. In this paper we have applied 
the discrete energy method to test problems commonly 
found in numerical relativity and analyzed the numeri- 
cal and long-term stability of the resulting schemes, with 
and without black hole excision. We have also investi- 
gated the numerical origin of fast growing modes which 
may artificially appear in numerical solutions, and used 
the energy method to discuss possible remedies. 

In particular, we have presented explicit examples of 
the techniques described in paper I and discussed what 
leads to ALTI and a way to remove or alleviate them. 
A message from this work is that it is not sufficient to 
construct numerically stable implementations, but also 
special care must be taken to devise suitable numerical 
methods to rule out, or alleviate ALTI. (For a somewhat 
related approach, see the recent work |3^). 

To illustrate how this is achieved we have adopted 
a series of examples, which although considerably sim- 
pler than Einstein equations, are relevant in the context 
of numerical relativity. In particular, in efforts to deal 
with black hole spacetimes and to understand the source 
of problems encountered with the gauge wave test. As 
shown throughout this work, even though there are ex- 
amples in which low-frequency ALTI are present and for 
which no sharp estimate is known, a careful analysis and 
grouping of the variables and/or the addition of (control- 
lable) artificial dissipation may in many cases be crucial 



APPENDIX A: CUBIC EXCISION FOR KERR 
BLACK HOLES 



Black hole excision removes the singularity by plac- 
ing an inner boundary on the computational domain. 
While the boundary removes the singularity, the ques- 
tion now becomes finding physically and mathematically 
consistent inner boundary conditions. Clearly, a general 
solution to this problem for an arbitrary boundary is un- 
known. The time-development of this boundary presents 
a further complication. Namely, if the boundary be- 
comes time-like, it will encounter the singularity in a fi- 
nite amount of time. These two issues can be resolved 
simultaneously by adopting an inner boundary with a 
space-like world tube. This ensures it will not hit the sin- 
gularity and that no boundary conditions are necessary 
if all the characteristics speeds are within the light cone, 
as the past domain of dependence of boundary points 
is wholly contained on the previous time slice. Some 
freedom remains in choosing the particular region to be 
excised, and clearly the largest possible excision bound- 
ary is preferable to avoid the strongest gradients of the 
spacetime. 

To illustrate some of these issues we now examine cu- 
bical excision for both the Schwarzschild and Kerr black 
holes, looking for the largest possible space-like excision 
cube. As we discuss below, the shape of the boundary 
plays a crucial role in satisfying the spacelike condition. 
We examine one face of the cube, a. x = const face. 
The condition that this surface is space-like implies that 
gxx ^ g^ ^g y^ jg ^^^ normal to the surface. Since the 
coordinates are fixed a-priori, we can only vary the posi- 
tion of the boundary to meet this condition. We find that 
this condition is quite difficult to satisfy, and possible for 
only very slowly spinning black holes. 

To see this in detail, consider the Kerr metric in Kerr- 
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Schild form |33- The quantity g^^ is 



9 



1 



ay 



where r is related to the Cartesian coordinates by 



(Al) 



Finding the minimum sized cube is more involved, as the 
faces of the cube need to be examined in addition to the 
corners l22l. 



y 



1 



(A2) 



The space-like character requirement on the surface im- 
plies that 



Mr^ 



r* + a'^z 



2~2 



rx -\- ay 

J.2 _|_ Q,2 



>1. 



(A3) 



For the Schwarzschild black hole, a — 0, r^ = x^- 
and the above condition becomes 

Mx^ 
2^^ > 1. 



-y^ + z^, 



(A4) 



The corner of the cube is the limiting point [y — z ~ x) 
for the largest cube, giving the bound x < 2Af/(3\/3), or 
X < 0.385M. 

For the spinning black hole, the analysis is slightly 
more involved. The limiting point for the largest cube 
is given hy y = —x, z = ±a;, giving the condition 



2MV2x3(V3x- 2-1/2^)2 4:^2 



22-1/2 



(3x2 



',2\2 



> 



(A5) 



Examining this condition shows that cubical excision can 
not be used when a > 0.082, as shown in Fig. [2] Fur- 
thermore, the singularity is not at a point, but is ring- 
like, implying that the cube may not be arbitrarily small. 



This limitation of cubical excision for the 
Schwarzschild spacetime was noted before |3^, by 
requiring that the excision boundary has no in-coming 
modes in a particular formulation of Einstein equations. 
A similar analyzis to that of |23 for the Kerr spacetime 
and cubical excision regions has also recently presented 
in [22|. Here we would like to point out the following. 
First, when the eigenmodes are complete and the speeds 
are physical (they lie within the null cone) those results 
and the one presented here are equivalent, as it is 
rooted in the physical nature of the excision boundary. 
Second, when dealing with a formulation admitting 
eigenspeeds outside the light cone, or a non-complete 
set of eigenmodes, the current straightforward analysis 
provides necessary conditions that have to be met. 
These necessary conditions are useful when dealing with 
any formulation, in particular one whose eigenspeeds 
might not be known. Finally, we emphasize that these 
results are coordinate dependent; nevertheless, for the 
families often used in numerical relativity tests, the same 
result applies. In particular for the family of slicings 
reported by Martel and Poisson |33 (which contains 
both the Kerr-Schild and Painleve-GuUstrand slicings 
of Schwarzschild) . 
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In order to obtain an estimate in the non-linear case, 
one might need to include derivatives of the fields in the 
energy as well. 

Despite the name, this energy does not need to corre- 
spond to the physical energy of the system. 
Though some energy estimate is still available, as dis- 
cussed in Paper I, and the scheme is at least numerically 
stable. 

The superscript index (1) is used to denote that this is 
the first energy defined for this system. 
This time denoted with the superscript (2). 
Note that the initial maximum amplitude of 11 is of order 
one despite the large factor in its definition, max{n} — 
1.25''/8 ^ 0.6. 



